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Fusion of very neutron rich nuclei may be important to determine the composition and heating of 
the crust of accreting neutron stars. We present an exploratory study of the effect of the neutron-gas 
environment on the structure of nuclei and the consequences for pycnonuclear fusion cross-sections in 
the neutron drip region. We studied the formation and properties of Oxygen and Calcium isotopes 
embedded in varying neutron-gas densities. We observe that the formed isotope is the drip-line 
nucleus for the given effective interaction. Increasing the neutron-gas density leads to the swelling 
of the nuclear density. We have used these densities to study the effect of this swelling on the 
fusion cross-sections using the Sao-Paulo potential. At high neutron-gas densities the cross-section 
is substantially increased but at lower densities the modification is minimal. 

PACS numbers: 25.60.Pj, 26.60.-l-Gj, 97.80.Jp 


I. INTRODUCTION 

Recent advances in radioactive beam technologies have 
opened up new experimental possibilities to study fu¬ 
sion of neutron rich nuclei [1]. Furthermore, near barrier 
fusion cross sections are relatively large so experiments 
are feasible with modest beam intensities. In addition, 
measurements are possible at the TRIUMF ISAC facility 
and in the near future at the NSCL ReA3-6 reacceler¬ 
ated beam facility. Other radioactive ion beam facilities 
include ATLAS-CARIBU at Argonne National Labora¬ 
tory, SPIRAL2 at GANIL (France), and RIBF at RIKEN 
(Japan). Note that the dynamics of the neutron rich skin 
of these nuclei can enhance the cross-section over that 
predicted by a simple static barrier penetration model. 
For example, neutrons may be transferred from the neu¬ 
tron rich beam to the stable target. Fusion of very neu¬ 
tron rich nuclei, near the drip line, raise very interesting 
nuclear structure and nuclear dynamics questions. 

Neutron stars, in binary systems, can accrete material 
from their companions. This material undergoes a vari¬ 
ety of nuclear reactions [2]. First at low densities, con¬ 
ventional thermonuclear fusion takes place, see for exam¬ 
ple [3]. Next at higher densities, the rising electron Fermi 
energy induces a series of electron captures [4] to produce 
increasingly neutron rich nuclei. Finally at high densities, 
these very neutron rich nuclei can fuse via pycnonuclear 
reactions. Pycnonuclear fusion is induced by quantum 
zero point motion [5, 6]. The energy released, and the 
densities at which these reactions occur, are important 
for determining the temperature and composition profile 
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of accreting neutron star crusts. The existence of the in¬ 
ner neutron-star crust, in which very neutron rich nuclei 
are immersed in a gas of neutrons raises the question, 
what is the impact of this neutron gas on nuclear fu¬ 
sion rates? This neutron drip region is believed to occur 
for densities in the approximate range of a few xlO^^ 
to 8 X 10^^ g/cm^. One of the early studies of the in¬ 
ner crust, consisting of drip-line nuclei combined with 
the background neutron gas, had been done by Negele 
and Vautherin [7]. Therefore understanding fusion re¬ 
actions of neutron rich isotopes near the drip line are 
important. Horowitz et al. [8] calculate the enhancement 
in fusion rates from strong ion screening using molecu¬ 
lar dynamics simulations, and find that can 

fuse near 10^^ g/cm^, just before neutron drip. Extensive 
studies of the astrophysical S{E) factors have been done 
using densities emanating from microscopic calculations 
and a barrier penetration model for fusion [9, 10]. Fur¬ 
thermore, this fusion can take place in the background 
neutron gas that is present in the inner crust of a neutron 
star. The possible effect of the neutron gas background 
was discussed in Ref. [10] by empirically changing the 
barrier height and width. Here, we study this effect by 
considering the presence of the background neutron gas 
microscopically by directly including the neutron gas and 
the nucleus in the same framework. We explicitly calcu¬ 
late the self-consistent proton and neutron densities of a 
single nucleus in equilibrium with the background neu¬ 
tron gas. Furthermore, for astrophysical applications, it 
seems clear that this adiabatic approach is the one that 
is relevant for calculating fusion rates in the inner crust. 
The neutron gas should have plenty of time to adjust to 
the presence of a nucleus. 

The paper is organized as follows. Our computational 
approach to consider the presence of neutron gas together 
with the nucleus and the model for calculating the fusion 
cross-sections is discussed in Sec. H. Computational re¬ 
sults for the Oxygen and Calcium systems is described 
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in Sec. III. Finally, these results are discussed and we 
conclude in Sec. IV. 


II. COMPUTATIONAL DETAILS 
A. Computational setup 


Hartree-Fock (HF) calculations were done in a three- 
dimensional Cartesian geometry with no symmetry as¬ 
sumptions and using the Skyrme effective nucleon- 
nucleon interaction [11]. The inhnite neutron star crust 
environment is simulated by using a three-dimensional 
Cartesian box with periodic boundary conditions for 
both the bound and neutron gas states as well as the 
solution of the Poisson equation for the Coulomb po¬ 
tential, which is performed using Fast-Fourier Transform 
(FFT) techniques [12]. The Coulomb solution assumes 
global neutrality which means that the proton charges 
are compensated by a homogeneous negative electron gas 
cloud. In practice, this is achieved by setting the zero- 
momentum part of the Coulomb field in Fourier space to 
zero. The code uses the basis-spline collocation method 
for the lattice discretization of the HF equations using pe¬ 
riodic boundary-conditions as described in Refs. [13-15]. 
The HF equations are solved using the damped gradi¬ 
ent iteration method. The Skyrme parametrization used 
was SLy4 [16]. In addition to providing a good descrip¬ 
tion of nuclei this interaction has been used to produce 
an equation of state for neutron stars [17]. 

For the choice of initial states to be used in HF min¬ 
imization we have tried a number of choices, which all 
resulted in the same identical solution. One can first gen¬ 
erate any isotope of the desired nucleus by solving the HF 
equations as described above and subsequently combine 
these states with a large number of free neutron gas states 
and minimize the entire system again. Alternately, one 
can simply choose a number of free proton states together 
with a large number of free neutron states and minimize 
this system. Both methods result in exactly the same nu¬ 
merical solution with a drip-line isotope corresponding to 
the nucleus with the given number of protons embedded 
in a given density of neutron gas states. Initial states, 
are spinors with a non-zero upper component in case 
of time-reversal invariance. In case of no time-reversal 
invariance the number of states are doubled by adding 
spinors having non-zero lower components as well. They 
satisfy the periodicity condition 

^n(l' + L) = V’nW , (1) 


where n = {nx,ny,nz) and Ua taking on integer values 
—Vq, ..., +Na. Free states satisfying the above period¬ 
icity condition are simple plane-wave states with the ap¬ 
propriate normalization 


V’„(r) = 


\/ L^LyL^ 


X + kri y + k„^ z) ^ 


( 2 ) 


where = 2'Kna/La and Xn is an up or down spinor. 
The initial neutron and proton densities are perfectly 
uniform filling the entire numerical box. These initial 
states comprise the total number of states used in the 
self-consistent HF problem using the Skyrme interaction. 
For even number of states time-reversal is valid and the 
HF single-particle Hamiltonian only depends on particle 
density, p, kinetic energy density, r, and the spin-orbit 
pseudotensor J through the single-particle states [16] 

= eA</>A A = 1,...,V. (3) 

As the HF iterations proceed (preserving orthogonality 
for the entire system) some states evolve to form a bound 
nuclear system while the others remain as gas states 
showing some non-uniformity due to the presence of shell 
effects. 


B. Fusion cross-sections 

The Sao Paulo model of fusion calculates an effective 
nuclear potential based on the density overlap between 
colliding nuclei [18, 19]. Sub-barrier fusion cross-sections 
can then be calculated via tunneling. The model can be 
easily applied to a very large range of fusion reactions 
and qualitatively reproduces many experimental cross- 
sections [20, 21]. Recently this model was used to tab¬ 
ulate astrophysical S factors describing fusion of many 
carbon, oxygen, neon and magnesium isotopes for use in 
astrophysical simulations [9], see also Ref. [22]. 

In this section we describe the Sao Paulo barrier pen¬ 
etration model to calculate fusion cross-sections. This 
starts with the double folding potential Vf{R) [18, 19], 

Vf{R)= J d^ri(fr2Piiri)p2{r2)Vo6{ri-r2-'R). (4) 

Here pi and p 2 are the densities of the two nuclei and 
Vo = —450 MeV-fm^. From Vf a nonlocal potential 
Vn{R,E) is constructed, Vn{R,E) = 
where v is the local relative velocity [18, 19] between the 
two nuclei at separation R (c is the speed of light) 

v^iR,E) = -[E-VciR)-VNiR,E)]. (5) 

Here p is the reduced mass and Vc{R) is the Coulomb 
potential at R. In practice, we use FFT techniques to 
calculate Vf{R) as well as the Coulomb potential Vc{R) 
(instead of using the point Coulomb formula). The ve¬ 
locity equation (5) has to be solved by iteration at each 
value of R and E. 

Note that the neutron gas background could behave 
differently for two nuclei in close proximity then it does 
for only a single nucleus. However for simplicity, in this 
first study, we consider only a single nucleus in the back¬ 
ground gas at a time in order to get the density profiles 
shown in Fig. 2. We then use these profiles in Eq. 4 
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and assume they are unmodified by the presence of the 
second nucleus. 

The fusion barrier penetrabilities are ob¬ 

tained by numerical integration of the two-body 
Schrodinger equation 




L{L + iW 

' 2/ri?2 


+ V{R,E)-E 


Z = o, 


( 6 ) 


using the incoming wave boundary condition (IWBC) 
method [23]. The potential V{R,E) is the sum of nu¬ 
clear and Coulomb potentials. IWBC assumes that once 
the minimum of the potential is reached fusion will occur. 
In practice, the Schrodinger equation is integrated from 
the potential minimum, where only an incoming 

wave is assumed, to a large asymptotic distance, where 
it is matched to incoming and outgoing Coulomb wave- 
functions. The barrier penetration factor, TL{Ec.m.) is 
the ratio of the incoming flux at i?min to the incoming 
Coulomb flux at large distance. Here, we implement the 
IWBC method exactly as it is formulated for the coupled- 
channel code CCFULL described in Ref. [23]. This gives 
us a consistent way for calculating cross-sections at above 
and below the barrier via 


af{E,,^,)^ + ■ (7) 

L=0 


III. RESULTS 

All the calculations presented here were done using a 
three-dimensional cubic Cartesian box with 31 fm sides 
and 1.0 fm lattice spacing. With the basis-spline method 
this gives highly accurate results for the HF problem [13]. 
We studied two systems with Z = 8 and Z = 20 protons 
inside a neutron gas. The resulting nuclei are always 
spherical and density profiles can be obtained by taking 
a cut along a particular axis. Our mesh includes the 
origin in all directions. We have repeated some of these 
calculations in a cubic bot with 41 fm sides and for the 
same neutron-gas density the results were numerically 
indistinguishable. 


A. Z = 8 system 

In these calculations we start with 8 proton states and 
a number of neutron states ranging from 50 — 1020. The 
Skyrme SLy4 force gives as the slightly bound drip¬ 
line nucleus in free-space. For the above range of neu¬ 
tron states we also find the to be the bound part 
of the system. As the number of neutron states is in¬ 
creased an overall negative potential is developed perme¬ 
ating the entire box. Figure 1 shows the neutron and 
proton mean-field potentials (solid lines) for Z = 8 and 
520 neutrons. The dashed lines indicate the energies of 
the bound single-particle states with degeneracies shown 



FIG. 1. (Color online) Mean-field potentials for neutrons and 
protons (solid lines) for the system with Z = 8 and 520 neu¬ 
trons. The dashed lines indicate the energies of the bound 
single-particle states with degeneracies shown in brackets. 


in brackets. We define Abound as the highest neutron s.p. 
state below the continuum threshold. Consequently, the 
bound and gas densites become 


.Abound 


Pbound — ^ ^ |0 a| : 

A=1 

(8) 

Pfree = ^ ^ |0 a| 

(9) 


.^^.^bound 


In this case the asymptotic value of the neutron po¬ 
tential is about —8 MeV. As it is the case in free-space 
increasing neutron number leads to deepening of the pro¬ 
ton potential. In Fig. 2 we show the density profiles for 
neutrons and protons as well as the total density as a 
function of the number of neutron-gas states. The top 
frame shows the total density behavior as the neutron- 
gas density is increased. The curves labeled n = 20 
correspond to free-space nucleus. As the external 
neutron-gas density is increased the bound system swells 
up in a way similar to a density scaling p{r) —)■ p{sr) 
with s < 1 as discussed in Ref. [24]. While the peak of 
the total density decreases from the free-space value of 
0.16 fm“^ to as low as 0.068 fm“^ for the 1000 external 
neutron-gas state case, the tail region flattens and de¬ 
velops a larger spatial extent, since the total integral re¬ 
mains to be 28. The density profiles are symmetric about 
a: = 0 and the numerical box extends to larger values then 
shown in the figure. Using these densities in Eq. (4) we 
have calculated the corresponding ion-ion folding poten¬ 
tials as well as the Coulomb potential. In Fig. 3 these 
potentials are plotted for a range of external neutron- 
gas densities. What is observed is that for neutron-gas 
densities in the range pgas = 2 — 4 x 10^^ gm/cm^ the 
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FIG. 2. (Color online) (a) Total density profiles for bound 
states; (b) density profiles for bound neutrons (solid lines) 
and protons (dashed lines), for the system with Z = 8 and 
n = 20, 520, and 1020 neutrons. 



FIG. 3. (Color online) Ion-Ion potentials Vsp{R) = Vf{R) + 
Vc{R) for isotope as a function of the external neutron- 
gas density. Also shown is the point-Coulomb interaction. 
Densities are in units of gm/cm^. 

effect of the gas is not changing the ion-ion potential in 
comparison to the free-space case in a considerable way. 
However, for gas densities above 10^^ gm/cm^ a very sig- 
nihcant change is observed. The free-space barrier has 
a peak value of 7.87 MeV located at i? = 10.7 fm. As 



FIG. 4. (Color online) Astrophysical S factor versus center of 
mass energy for fusion of isotope as a function of external 
neutron-gas density. Cross-sections are calculated using the 
Sao Paulo barriers and the IWBC method. 

the external gas density is increased the corresponding 
barrier height is reduced to 7.69, 7.23, and 6.39 MeV 
with the peak location moving outward at 10.8, 11.7, 
and 13.3 fm. This behavior is very similar to what is 
observed in free-space as one goes up in the oxygen iso¬ 
tope chain [25]. Unfortunately, the dynamical density 
constrained time-dependent Hartree-Fock (DC-TDHF) 
methodi [26-28] used in Ref. [25] is not applicable in this 
situation since it requires a fully dynamical calculation. 

These together with the calculation of the energy de¬ 
pendence from Eq. (5) allows the calculation of fusion 
cross-sections as a function of external neutron-gas den¬ 
sity, using the IWBC method discussed in the previous 
section. Figure 4 demonstrates the effect of neutron-gas 
seen in the potentials on the astrophysical S factor. The 
S-factor for different external neutron-gas densities start 
to deviate from each other as soon as the center-of-mass 
energy falls below the barrier. Even for the lowest gas 
density of 2.9 x 10^^ gm/cm^ the difference with the free- 
space value is about a factor of two at the center-of-mass 
energy of 2 MeV. The difference at higher gas densities 
are about 1 — 3 orders of magnitude larger then the free- 
space values at sub-barrier energies. 

In Fig. 5 we plot the cross-sectional density (j/ = 0 
plane) profile of the ^®O-|-500n system. In general we see 
a higher gas density in the vicinity of the nucleus as can 
also be seen as the dotted line in Fig. 2(a). The energy 
of this low density neutron gas is so low that even very 
small ’’shell effects” can lead to nonuniform densities. In 
some of the cases we studied the neutron density is low 
enough to have shell structures leading to the formation 
of neutron arms seen in Fig. 5. If we replicate this cubic 
box in three-dimensions one sees a lattice like structure 
linked by these neutron arms. These effects are driven 
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FIG. 5. (Color online) The cross-sectional density (i/ = 0 
plane) profile of the ^*O-|-500n system. The neutron density 
is low enough to have shell structures leading to the formation 
of neutron arms. 


FIG. 6. (Color online) Mean-field potentials for neutrons and 
protons (solid lines) for the system with Z — 20 and 540 
neutrons. The dashed lines show the mean-field potential of 
®°Ca in free-space. 


mainly by the periodic boundary conditions here. Irreg¬ 
ularities in the grid of nuclei will produce more irregular 
structures of these arms. Since the density in the arms is 
of the order of free density such that they do not affect 
the analysis of the bound part. 

B. Z = 20 system 

We have repeated the same study by starting with 20 
proton states and a number of neutron states ranging 
from 140 — 1040. Using the Skyrme SLy4 we get ®°Ca 
as the slightly bound drip-line nucleus in free-space. For 
the above range of neutron states we also find the ®°Ca 
to be the bound part of the system. As the number 
of neutron states is increased an overall negative poten¬ 
tial is developed permeating the entire box. Figure 6 
shows the neutron and proton mean-field potentials (solid 
lines) for Z = 20 and 540 neutrons. The dashed lines 
show the mean-field potentials in free-space. In this case 
the asymptotic value of the neutron potential is about 
—6.5 MeV. At higher neutron densities we do observe 
the tendency for ^^Ca to be the drip-line nucleus but the 
tendency is very weak and for practical purposes consid¬ 
ering ®°Ca is sufficient for the general purposes of this 
study. 

In Fig. 7 we plot the density profiles for neutrons and 
protons as well as the total density as a function of the 
number of neutron-gas states. The top frame shows the 
total density behavior as the neutron-gas density is in¬ 
creased. The curves labeled n = 40 correspond to free- 
space ®°Ca nucleus. As the external neutron-gas den¬ 
sity is increased the bound system swells up as in the 


Z = 8 case. While the peak of the total density de¬ 
creases from the free-space value of 0.165 fm“^ to as low 
as 0.048 fm“^ for the 1000 external neutron-gas state 
case, the tail region flattens and develops a larger spatial 
extent, since the total integral remains to be 60. The 
density profiles are symmetric about x = 0 and the nu¬ 
merical box extends to larger values then shown in the 
figure. Corresponding ion-ion folding potentials calcu¬ 
lated by using these densities in Eq. (4) are plotted in 
Fig. 8 for a range of external neutron-gas densities. Again 
what we observe is that for neutron-gas densities in the 
range pgas = 2 — 4x 10^^ gm/cm^ the effect of the gas in 
changing the ion-ion potential compared to the free-space 
case is not significant. However, for gas densities above 
10^^ gm/cm^ a very significant change is observed. The 
free-space barrier has a peak value of 47.4 MeV located 
at i? = 11.3 fm. As the external gas density is increased 
the corresponding barrier height is reduced to 46.3, 42.9, 
and 37.4 MeV with the peak location moving outward at 
11.6, 12.5, and 14.3 fm. Figure 9 shows the fusion cross- 
sections calculated using the ion-ion potentials shown in 
Fig. 8. The dramatic rise of the cross-section is obvious 
as the neutron gas density becomes significantly higher 
than the minimum neutron drip density. 

IV. SUMMARY AND DISCUSSION 

Pycnonuclear reactions are expected to occur between 
very neutron-rich nuclei and in a dense background of of 
a neutron gas [10]. We have made an exploratory study 
of the effect of the external neutron gas on nuclear fusion 
in the regime between the start of the neutron drip region 
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FIG. 7. (Color online) (a) Total density profiles for bound 
states; (b) density profiles for bound neutrons (solid lines) 
and protons (dashed lines), for the system with Z = 20 and 
n = 40, 540, and 1040 neutrons. 


ment of the Coulomb potential with periodic boundary 
conditions. For the calculation of fusion cross-sections 



FIG. 9. (Color online) Fusion cross-section as a function 
of center-of-mass energy for the fusion of two ®'^Ca isotopes 
calculated for changing external neutron-gas density. Cross- 
sections are calculated using the Sao Paulo barriers and the 
IWBC method. 



FIG. 8. (Color online) Ion-Ion potentials Vsp{R) = Vf{R) + 
Vc{R) between two ®°Ca isotopes as a function of the exter¬ 
nal neutron-gas density. Also shown is the point-Coulomb 
interaction. Densities are in units of gm/cm®. 

and the melting region. For our study we have used an 
approach that treats the nuclei and the extra neutrons in 
a unified manner. The computations are done in full 3D 
with periodic boundary conditions and a special treat- 


we have used the Sao-Paulo model. In our calculations 
we observe that for lower background densities the cross- 
sections do not change in a very significant manner. On 
the other hand as we increase the neutron gas density 
we observe the swelling of the nuclei that results in the 
lowering of the ion-ion potential barriers and significant 
increase in the fusion cross-sections. At densities higher 
than the neutron-drip regime the melting of the nuclei 
can be observed. 

While our methods give us a good understanding of fu¬ 
sion under these conditions more precise computations, 
including the full effects of pairing and effective inter¬ 
actions tailored for neutron star crust [29], may reduce 
some of the observed shell effects and modify some of the 
results but the main features observed are expected to 
remain the same. In addition, movement of nuclei inside 
the neutron gas as they approach each other may cause 
ripples and waves in the neutron background that can 
also influence these results. However, in the adiabatic 
limit these effects are not expected to be very large. 
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